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Propellant-free  Spacecraft  Relative  Maneuvering  via  Atmospheric  Differential  Drag 


FINAL  PERFORMANCE  REPORT 

Abstract 

At  low  Earth  orbits  (LEO),  a  differential  in  the  drag  acceleration  between  spacecraft  can  be  used  to 
control  their  relative  motion.  This  differential  allows  for  a  propellant -free  method  for  performing  relative 
maneuvers,  which  can  reduce  the  costs  of  spacecraft  formations.  In  this  project  atmospheric  differential 
drag  (DD)  based  nonlinear  controllers  are  presented  that  can  be  used  for  planar  relative  maneuvers  of  two 
spacecraft.  The  controllers  were  tested  using  high  fidelity  Systems  Tool  Kit  (STK)  simulations  for  re¬ 
phase,  fly-around,  and  rendezvous  maneuvers.  Furthermore,  the  atmospheric  density  varies  in  time  and  in 
space  as  the  spacecraft  travel  along  their  orbits.  In  this  project  a  localized  density  predictor  based  on 
Neural  Networks  (NNs)  was  also  developed.  The  predictor  uses  density  measurements  or  estimates  along 
the  past  orbits  and  can  use  a  set  of  proxies  for  solar  and  geomagnetic  activities  to  predict  the  value  of  the 
density  along  the  future  orbits  of  the  spacecraft.  The  performance  of  the  localized  predictor  is  studied  for 
different  NN  structures,  testing  periods  of  high  and  low  solar  and  geomagnetic  activity  and  different 
prediction  windows.  Comparison  with  previously  developed  methods  show  substantial  benefits  in  using 
neural  networks,  both  in  prediction  accuracy  and  in  the  potential  for  spacecraft  onboard  implementation. 
List  of  symbols 

a,  b,  c,  d  Constants  in  Schweighart  and  Sedwick  model 

aDrei  Magnitude  of  the  aerodynamic  drag  acceleration  experienced 

by  spacecraft 

aDcrit  Magnitude  of  the  differential  drag  acceleration  than  ensures 

stability 

PreI  Differential  accelerations  caused  by  any  additional  perturbations 

Abc  Element  of  matrix  Aj  to  which  aDcrit  is  the  most  sensitive 

Aj  Stable  linear  reference  state  space  matrix 

A,  Schweighart  and  Sedwick  model  state  space  matrix 

B  Spacecraft’s  ballistic  coefficient,  also  bias  in  the  context  of  the 

neural  network  predictors 
C[>  Spacecraft  drag  coefficient 

D  Number  of  delays  in  the  hidden  layer 

Dst  Geomagnetic  index 

e  Tracking  error  vector 
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Nonlinearities  in  the  spacecraft  dynamics,  or  nonlinearities  in  general  in  the 
context  of  Lyapunov  stability,  also  nonlinear  function  in  the  context  of  the  neural 
network  predictors 
Solar  index 

Target’s  initial  orbit  inclination 

Second-order  harmonic  of  Earth  gravitational  potential  field  (Earth 
flattening) 

Lower  and  upper  bounds  for  the  optimization 
Spacecraft’s  mass 
Mean  motion 

Number  of  samples  in  a  data  set 

Solution  matrix  of  the  Lyapunov  equation  and  its  vectorized  form 
Selected  Lyapunov  equation  matrix  and  its  vectorized  form 
Pearson  correlation  coefficient 
Earth  mean  radius 

Spacecraft  cross-wind  section  area  for  chaser  and  target  spacecraft 

Sampling  period  of  the  density  data 

Matrix  transformations 

Feedback  linear  control  for  reference  system 

On-off  control  signal 

Permutation  matrices 

Lyapunov  function 

Prediction  window 

State  space  vector  including  relative  position  and  velocity  of  the 
spacecraft  system  in  the  LVLE1  orbital  frame 
Linear  reference  model  state  space  vector  in  the  LVLH  orbital 
frame 

Ballistic  coefficient  differential 

Magnitude  of  the  modification  made  to  matrix  Ad using  the 
switching  adaptation 

Modifications  made  to  matrix  At!  for  the  optimized  adaptation 
Neural  network  input 
Gravitational  parameter 
Neural  network  output 


Atmospheric  density 
Predicted  Atmospheric  density 
Neural  network  target  value 


1.  Introduction 


Leonard  et  al.  [1]  introduced  the  concept  of  DD  as  an  alternative  method  for  generating  the  control 
forces  required  by  the  relative  maneuvers  on  the  orbital  plane  at  LEO.  The  method  consists  of  varying  the 
aerodynamic  drag  experienced  by  different  spacecraft,  thus  generating  differential  accelerations  between 
them.  The  interest  towards  this  methodology  comes  from  the  desire  for  efficient,  autonomous  spacecraft 
relative  maneuvering.  To  increase  the  efficiency  and  economic  viability  of  such  maneuvers,  propellant 
consumption  must  be  reduced.  Furthermore,  since  there  is  no  propellant  exhaust,  sensitive  onboard 
sensors  will  not  be  affected  during  maneuvers.  However,  maneuvering  using  DD  requires  the  ability  to 
generate  magnitudes  of  differential  acceleration  larger  than  the  accelerations  caused  by  any  other 
perturbations;  hence,  the  spacecraft  must  be  able  to  change  their  ballistic  coefficients  by  the  necessary 
amounts  and  their  orbits  must  be  low  enough  for  providing  enough  drag  force.  Additionally,  DD  results  in 
increased  orbital  decay. 

Using  DD  for  spacecraft  relative  maneuvering  presents  two  main  challenges.  Firstly,  motion  on  the 
orbital  plane  is  to  be  controlled;  therefore,  two  degrees  of  freedom,  motion  along  the  orbit  and  along  the 
radial  direction  of  the  spacecraft,  are  to  be  controlled.  Conversely,  drag  force  acts  only  along  the  orbit  of 
the  spacecraft,  and  influences  the  motion  along  the  radial  direction  through  the  coupling  in  the  dynamics. 
Accordingly,  the  spacecraft  system  is  under  actuated  when  using  DD  for  controlling  its  dynamics  on  the 
orbital  plane.  Secondly,  the  magnitude  of  the  DD  acceleration  fluctuates  during  the  maneuver  as  the 
spacecraft  encounters  regions  of  the  thermosphere  with  varying  atmospheric  density.  These  variations  are 
difficult  to  predict  on  board.  This  means  that  the  control  force  (drag  force)  is  difficult  to  know  with 
certainty  prior  to  the  maneuver.  This  works  presents  methods  specially  developed  for  addressing  these 
issues. 

This  work  presents  the  development  of  novel  nonlinear  control  strategies,  based  on  the  Lyapunov 
approach,  to  perform  spacecraft  relative  maneuvers  at  LEO,  exploiting  atmospheric  drag  forces. 
Furthermore,  a  localized  predictor  based  on  NNs  for  the  estimation  of  future  values  of  the  drag 
acceleration  experienced  by  spacecraft  along  its  orbit  is  also  developed  as  complement  for  the  Lyapunov 
controllers.  The  use  of  the  predictor  will  provide  the  guidance  and  control  systems  with  estimates  of  the 
available  drag  force  along  the  future  orbit  of  the  spacecraft.  Such  tools  will  provide  a  framework  from 
which  autonomous  propellant-less  spacecraft  relative  maneuvering  using  atmospheric  DD  acceleration 
can  be  realized. 

The  problem  of  designing  a  control  system  for  the  maneuvering  using  DD  becomes  the  problem  of 
designing  a  real-time  logic  (i.e.,  the  control  action  is  computed  as  the  maneuver  progresses)  to  command 
the  deployment  or  retraction  of  the  drag  surfaces,  with  the  intent  of  forcing  the  satellites  to  follow  a 
desired  trajectory,  a  reference  model  or  simply  regulate  to  a  desired  final  state.  The  sought-for  logic  needs 


to  be  based  on  the  assumption  that  the  control  is  either  positive  maximum,  negative  maximum,  or  zero, 
neglecting  the  time  required  to  deploy  or  retract  the  drag  surfaces.  In  essence,  a  Lyapunov  function  of  the 
tracking  error  is  selected,  and  the  control  signal  is  chosen  so  that  the  tracking  error  converges  to  zero  (i.e., 
the  first  order  time  derivative  of  the  Lyapunov  function  is  negative),  thus,  the  nonlinear  dynamics  of  the 
system  (that  is  the  real  motion  of  the  spacecraft)  are  forced  to  follow  a  desired  trajectory,  reference  model 
or  desired  final  state  (regulation).  This  simplifies  the  control  problem,  since  the  desired  trajectory  can  be 
designed  using  controlled  linear  dynamics  approximating  the  reality  of  spacecraft  relative  motion.  The 
controllers  developed  by  the  author  build  upon  and  improve  work  presented  and  tested  in  [2-4], 

In  order  to  enhance  the  performance  of  the  Lyapunov  controller,  a  way  for  adapting  the  Lyapunov 
function,  in  terms  of  the  amount  of  drag  acceleration  necessary  for  stability,  was  developed.  The 
definition  of  appropriate  Lyapunov  functions  is  a  challenge  that  varies  from  problem  to  problem,  and  a 
widely  studied  theory  exists  (see  [5-7]).  In  this  work,  a  quadratic  Lyapunov  function  of  the  tracking  error 
is  defined,  and  its  positive  definite  matrix  is  changed  using  adaptations,  effectively  changing  the 
Lyapunov  function  in  real-time  to  improve  controller  performance  (control  effort  and  maneuver  duration) 
during  the  maneuvers.  Two  adaptation  strategies  were  developed,  the  switching  adaptation  and  the 
optimized  adaptation.  The  adaptations  are  achieved  through  analytical  expressions  giving  the  dependence 
of  the  critical  value  (minimum  magnitude  of  the  DD  acceleration  necessary  to  retain  Lyapunov  stability) 
on  the  chosen  stable  linear  model.  By  means  of  these  relationships,  the  amount  of  DD  acceleration 
necessary  to  maintain  stability  is  reduced  while  maneuvering,  achieving  a  better  control  authority  margin 
in  real-time.  Overall,  the  development  of  the  adaptive  Lyapunov  controllers  provides  a  valuable  method 
that  can  be  implemented  onboard  real  spacecraft,  even  small  spacecraft,  with  limited  computing 
capabilities.  In  fact,  the  Lyapunov  controllers  require  only  onboard  measurements  (relative  position  and 
velocity)  that  would  be  available  during  flight. 

To  develop  the  density  predictor,  a  similar  approach  to  that  of  Stastny  et  al.  [8]  was  used  in  this  work. 
However,  NNs  are  used  instead  of  a  linear  model.  NNs  are  capable  of  forecasting  nonlinear  behaviors 
since  they  contains  nonlinearities  in  their  neurons,  and  therefore  it  have  the  potential  to  accurately  model 
the  nonlinear  behavior  of  the  density  along  the  orbit  of  the  spacecraft.  To  train,  validate,  and  test  the  NNs, 
density  data  from  the  CHAllenging  Minisatellite  Payload  (CHAMP,  see  [9])  mission  was  used. 

The  main  contributions  of  this  projects  to  the  state  of  the  art  are: 

1)  Analytical  expressions  for: 

a.  A  DD  control  law,  which  dictates  the  activation  of  the  drag  surfaces,  based  on  Lyapunov 
principles  (presented  in  [10]). 

b.  The  magnitude  of  the  DD  acceleration  that  ensures  stability  in  the  sense  of  Lyapunov  for 
the  system  (critical  value),  presented  in  [11,  12].  Paper  [11]  was  the  winner  of  the  best 


student  paper  award  for  the  category  Spacecraft  GNC  at  the  1  st  International  Academy  of 
Astronautics  Conference  on  Dynamics  and  Control  of  Space  Systems  and  was  later 
published  in  the  Acta  Astronautica  Journal  as  [12]. 
c.  The  partial  derivatives  of  the  critical  value  of  the  DD  acceleration  in  terms  of  Q 
(Lyapunov  equation  matrix),  and  4/  (reference  linear  dynamics  matrix).  These 
derivatives  were  first  introduced  for  the  case  of  regulation  only  in  [11,  12],  and  later 
generalized  for  tracking  a  trajectory  or  a  reference  model  in  [13,  14].  Matrices  Q  and  A_,i 
are  independent  variables  that  are  chosen  by  the  designer  in  the  development  of  the 
Lyapunov  control  law. 

2)  Adaptation  strategy,  based  on  the  partial  derivative  of  the  critical  value  in  terms  of  matrix  • 
which  chooses  in  real-time  an  appropriate  positive  definite  matrix  P  in  a  quadratic  Lyapunov 
function,  such  that  it  reduces  the  critical  value  (presented  by  the  author  in  [11,  12]). 

3)  A  second  adaptation  strategy  also  based  on  the  same  partial  derivative,  which  chooses  an 
appropriate  positive  definite  matrix  P  in  a  quadratic  Lyapunov  function,  such  that  it 
minimizes  the  critical  value. 

4)  Demonstration  of  feasibility  of  the  approach  via  numerical  high  fidelity  orbital  simulations 
using  STK  for  three  different  maneuvers:  a  re -phase,  a  fly-around  and  a  rendezvous 

5)  Development  of  neural  network-based  localized  models,  that  are  capable  of  forecasting  the 
density  to  be  encountered  by  a  spacecraft  along  its  orbit  for  prediction  windows  of  one,  eight 
and  32  orbits  into  the  future  (i.e.  approximately  90min,  12hr  and  two  days  respectively), 
presented  by  the  author  in  [15]. 

6)  Appropriate  design  of  the  NN  structure  using  different  parameters  such  as  the  sampling  rate 
of  the  data,  the  number  of  neurons  in  the  hidden  layer  and  the  number  of  delays  of  the  input. 

7)  Tests  of  the  NN  predictors  over  periods  of  high  and  low  solar  and  geomagnetic  activities. 

8)  Comparison  of  the  results  of  the  NN  predictors  with  a  simple  persistence  model,  a  linear 
model,  JB2006,  and  HASDM  (the  latter  three  obtained  from  [8])  for  the  one-orbit  forecast. 

The  following  assumptions  were  made  in  the  developments  of  this  project: 

•  Two  spacecraft  (target  and  chaser)  are  considered,  which  have  the  ability  of  changing  their 
ballistic  coefficients  by  deploying  or  retracting  a  set  of  surfaces. 

•  The  discussion  in  this  work  will  be  limited  to  in-plane  motion,  assuming  that  no  out-of-plane  is 
present,  or  that  it  is  controlled  with  different  means. 

•  The  drag  surface  deployment/retraction  time  is  assumed  to  be  negligible  with  respect  to  the 
duration  of  the  maneuvers.  Thus,  it  is  assumed  the  drag  surfaces  deploy/retract  instantly, 
generating  a  bang-off-bang  control  profile,  as  suggested  in  [2,  3,  16  and  17,]. 


•  The  attitude  is  stabilized  by  other  means  than  DD. 

This  document  is  organized  as  follows.  Section  2  presents  all  the  original  developments  in  this  project, 
including:  the  Lyapunov  controller,  the  adaptive  and  adaptive  optimized  Lyapunov  controllers,  and  the 
NN  based  predictors.  Section  3  summarizes  and  comments  the  results  of  this  work,  including:  the  results 
from  the  STK  simulations  of  the  relative  maneuvers  using  the  controllers  developed  and  the  results  of  the 
tests  done  to  the  NN  predictors.  Finally,  section  4  draws  the  conclusions. 


2.  Procedure 


In  this  work  the  spacecraft  relative  motion  is  represented  in  the  local  vertical  local  horizontal  (LVLH) 
reference  frame.  In  this  frame  the  origin  of  this  frame  orbits  with  the  target  spacecraft,  x  points  from  Earth 
to  the  tai'gct  spacecraft  (virtual  or  real),  z  points  along  the  angular  momentum  vector  of  the  target’s  orbit, 
and  y  completes  the  right-handed  frame.  The  following  equations  were  used  for  represent  the  relative 
motion  dynamics  for  a  chaser  and  target  spacecraft  in  the  LVLH  frame: 
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The  magnitude  of  the  relative  acceleration  caused  by  the  differential  aerodynamic  drag  for  the 
spacecraft  system  (target  and  chaser)  is  given  as: 


aDrel  ~ 


(2) 


where  AB  is  the  difference  in  ballistic  coefficients  between  the  target  and  chaser.  The  ballistic  coefficient 
is  defined  as 
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Furthermore,  the  Schweighart  and  Sedwick  [18]  linear  model  was  used  to  create  the  reference  model 
for  the  development  of  the  controller.  This  model  assumes  that  the  target’s  orbit  is  circular,  that  the  only 
perturbation  acceleration  is  caused  by  the  J2  perturbation,  and  that  the  separation  between  spacecraft  is 
small  compared  to  the  radii  of  their  orbits.  Schweighart  and  Sedwick  model  can  be  formulated  as  the 
following  system: 
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In  essence  the  Lyapunov  controller  used  in  this  work  forces  a  nonlinear  model  to  follow  a  stable  linear 
reference  system.  A  quadratic  function  of  the  tracking  error  (difference  between  the  state  of  the  nonlinear 
and  the  stable  linear  reference  systems)  is  defined.  Afterwards  a  control  law  aimed  to  satisfy  the 
Lyapunov  theorem  as  defined  in  [19].  The  stable  linear  reference  system  can  be  tracking  a  desired 
guidance  trajectory  or  can  be  regulated;  furthermore  the  stable  linear  reference  system  can  be  replaced  by 
a  desired  guidance  trajectory  or  by  a  constant  state  vector. 

V  =  eTPe,  e  =  xn-xd, 

(5) 

V  =  -erge  +  2 (eT PB a Drelu  -  erP(Arfx  -  /  (x)  +  Bud )) 

where  P  is  a  symmetric  positive  definite  matrix  and  e  is  the  tracking  error  vector,  u  is  the  command 
sent  to  the  surface  actuators,  the  matrix  Q_  is  chosen  such  that  a  Lyapunov  equation  is  satisfied 
(A/P+PA^-Q),  and  the  matrices  A d  and  B  represent  a  stable  linear  reference  dynamics  (in  this  work  the 
Schweighart  and  Sedwick  stabilized  using  an  LQR). 

The  resulting  control  law  presented  in  [10]  can  be  expressed  as: 

u  =  - signie 7  PB)  (6) 


It  is  worth  emphasizing  that  all  the  components  in  this  control  law  would  be  available  in  real-time 
onboard  a  spacecraft.  For  the  actual  implementation  of  the  controller  this  activation  strategy  is  applied 
every  10  minutes,  to  give  the  drag  forces  enough  time  to  change  the  net  accelerations  of  the  spacecraft. 

A  critical  value  (, aDcrit )  of  differential  drag  that  is  needed  to  maintain  stable  Lyapunov  control  was 
developed  in  [1 1,  12].  This  critical  value  is  given  as: 
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The  analytical  expression  for  the  partial  derivative  of  the  critical  value  with  respect  to  the  matrix  Aj  (in 
Equation  (8))  was  developed  in  [11,  12].  The  matrix  Aj  is  set  by  design,  so  it  can  be  manipulated  to 
reduce  the  value  of  the  critical  value  during  the  maneuver,  provided  that  it  remains  Hurwitz.  The  element 
of  matrix  Aj,  to  which  aDcri,  is  the  most  sensitive  (Abc,  the  entry  with  the  largest  partial  derivative)  is 
identified  by  calculating  the  partial  derivative. 


For  the  details  on  the  derivation  of  Equation  (8)  refer  to  [1 1,  12]. 


2.1  Adaptive  and  optimized  adaptive  Lyapunov  controllers 


To  improve  the  performance  of  the  Lyapunov  controller  an  adaptation  and  an  optimized  adaptation  of 
the  drag  surface  activation  strategy  (Equation  (6))  were  developed.  Both  adaptations  rely  (aDcrit)  and  the 
partial  derivative  (Equations  (7)  and  (8)). 


2.1.1  Adaptation  strategies 


Changing  the  elements  of  matrix  Aa  will  change  the  critical  value,  which  determines  how  much 
magnitude  of  the  relative  drag  acceleration  is  needed  to  attain  Lyapunov.  Furthermore,  modifying  this 
matrix  will  also  affect  the  control  law  via  the  matrix  P  (see  Equation  (6))  thus  changing  how  the 
Lyapunov  controller  reacts  to  the  tracking  error.  This  suggest  that  by  changing  this  matrix  the 
performance  of  the  controller  can  be  influenced.  Two  different  adaptations  were  developed  to  enhance  the 
performance  of  the  Lyapunov  controller  taking  advantage  of  this  feature.  These  adaptations  seek  to 
reduce  the  critical  value  via  modifications  to  matrix  A d  as  the  maneuver  progresses.  Both  adaptations  are 
implemented  as  illustrated  in  Figure  1 .  The  adaptations  are  implemented  such  that  the  modified  Aj  matrix 
is  still  stable  and  are  applied  at  the  same  time  as  the  activation  strategy,  which  is  every  10  minutes. 
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Figure  1.  Control  Strategy  diagram,  the  text  in  bold  represents  the  onboard  control  calculations,  occurring  every  10  minutes. 

By  calculating  the  partial  derivative  (Equation  (8))  the  entry  of  the  matrix  Ad.  to  which  aDcrit  is  the 
most  sensitive,  is  identified  (that  entry  which  has  the  partial  derivative  with  the  largest  magnitude),  which 
is  defined  as: 
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Once  this  entry  is  identified,  it  is  switched  to  a  small  value  (SA  =10 6  which  is  in  the  same  order  of 
magnitude  of  the  other  entries  of  A,).  The  sign  of  this  modification  is  chosen  such  that  it  reduces  the 
derivative  of  aDcri„  thus  inducing  a  downward  trend  in  the  behavior  of  the  critical  value  for  the  magnitude 
of  the  differential  acceleration.  By  reducing  this  critical  value,  the  overall  robustness  of  the  controller  is 
improved  since  the  control  margin  (the  difference  between  the  actual  value  of  the  DD  acceleration  and  the 
critical  value)  is  increased.  The  adaptive  variations  in  Aj  are  expressed  as: 
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where  ka  is  defined  by: 
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The  controller  resulting  from  applying  this  adaptation  to  the  Lyapunov  control  law  developed  in  the 
previous  section  is  called  adaptive  Lyapunov  controller. 

A  further  refinement  of  the  adaptation  strategy  was  achieved  by  including  an  optimization  method, 
more  specifically  MATLAB’s  fmincon  function,  in  the  adaptation.  The  f 'mine on  tool  is  an  optimization 
method  that  finds  the  minimum  of  a  constrained  multivariate  function. 

For  this  adaptation  the  partial  derivative  defined  in  is  again  calculated  and  just  as  in  the  switching 
adaptation,  the  element  of  matrix  to  which  a i)crit  is  the  most  sensitive  (Ahr)  is  found.  Afterwards,  the 
fmincon  function  is  used  to  minimize  aDcrit  in  terms  of  Ahc.  The  optimization  problem  solved  by  fmincon 
can  be  formulated  as: 


minimize  aDcnt(Alx) 

Abc 

subject  to  lb  <  Abc  <  up 


(12) 


where  and  lb  and  up  are  the  lower  and  upper  bounds  for  Ahc.  These  bounds  were  chosen  to  be  -10 6  and 
10 6  which  are  in  the  same  order  of  magnitude  of  the  other  entries  of  Aj.  The  solution  of  the  optimization 
problem  gives  the  sign  and  magnitude  (SAOp)  with  which  Abc  is  modified.  The  adaptive  variations  in  Ad  are 
expressed  as: 
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where  ka  and  is  defined  in  Equation  (11) 

It  should  be  noted  that  the  adaptations  occur  every  10  minutes  and  that  that  for  a  bang-off-bang  control 
the  At  from  Equations  (10)  and  (13)  is  essentially  zero. 


2.2  Density  Predictors 


The  NN  density  predictor,  developed  in  this  work,  uses  density  measurements  or  estimates  on  a  given 
orbit  and  a  set  of  proxies  for  solar  and  geomagnetic  activities  to  predict  the  value  of  the  density  along  the 
future  orbit  of  the  spacecraft. 


2.2.1  Neural  network  used 


A  time-delay  feed-forward  Neural  Network  structure  was  chosen  for  the  development  of  the  density 
predictor.  This  NN  architecture  does  not  include  any  feedback  loops,  hence  the  feed-forward  part  of  its 
name.  This  NN  architecture  contains  a  set  of  delays  at  the  input  layer  that  allow  for  retention  of  the 
evolution  of  the  inputs  in  time,  and  enhances  the  ability  of  the  network  for  forecasting  applications. 
Furthermore,  the  NN  predictor  contains  two  layers  (hidden  or  input  layer,  and  output  layer).  The  output 
layer  contains  one  single  linear  neuron.  The  number  of  neurons  and  delays  in  the  hidden  layer  were 
determined  by  testing  different  configurations.  The  results  of  these  experiments  are  included  in  the  results 
shown  in  section  3.  A  diagram  of  a  time-delay  feed-forward  Neural  Network  is  shown  in  Figure  2. 


Figure  2.  Diagram  for  a  FTDNN  with  two  layers,  two  inputs  (one  with  two  delays  and  the  second  one  with  one),  two  nonlinear 
neurons  in  the  input  layer  and  one  linear  neuron  in  the  output  layer. 


In  this  application,  the  input  to  the  NN  is  the  present  value  of  the  density  and  the  output  is  the 
predicted  value  over  a  predefined  prediction  window.  Additional  inputs  such  as  the  current  values  of  the 
solar  and  geomagnetic  indices  (Dst  and  F10.7)  can  also  be  included.  The  inputs  are  delayed  a  defined 
number  of  times  inside  the  NN  in  order  to  capture  some  of  their  time  evolution.  Such  formulation  is 
shown  in  the  following  expression: 
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where  /  is  the  overall  nonlinear  function  of  the  neural  network,  p  is  the  measured  density,  p  is  the 
predicted  density  value  (from  NN  output),  Wp  is  the  prediction  window,  ts  is  the  sampling  period  of  the 
data,  t  is  the  time,  and  I)  is  the  number  of  delays  in  the  hidden  layer. 

The  Levenberg-Marquardt  algorithm  (see  [20,  21])  was  used  to  train  the  NNs.  This  algorithm,  which  is 
included  in  MATLAB's  Neural  Network  Toolbox  ([22]),  was  chosen  since  it  often  has  higher  rates  of 
convergence  than  the  other  algorithms  provided  in  the  Toolbox.  The  mean  squared  error  (MSE),  as 
explained  in  Equation  (15),  was  selected  as  the  performance  function. 

MS£  =  -X(0,-r,)2  (15) 

ns  i= 1 

where  ns  is  the  number  of  samples,  and  T  is  the  target  value 

2.2.2  Density  data  used 

The  use  of  NNs  requires  data  sets  for  training,  validation,  and  testing  the  model’s  perfomiance.  The 
training  and  validation  sets  must  contain  data  covering  the  different  behaviors  to  be  modeled  by  the  neural 
network.  CHAMP  density  data  was  used  for  this  work.  Measurements  from  the  onboard  accelerometer 
onboard  the  CHAMP  satellite  allow  for  the  estimation  of  the  density  data.  These  data  is  available  online 
and  was  obtained  from  [23]. 

For  each  NN  the  training  density  data  was  divided  into  two  segments:  one  segment  of  past  values, 
assumed  to  be  available  for  the  training  and  validation  of  the  neural  network;  and  one  segment  of  future 
values,  which  are  values  of  density  that  would  not  be  available  during  training  and  validation,  but  instead 
are  used  exclusively  for  testing  the  NNs.  The  past  values  were  sampled  randomly  and  70%  were  used  for 
training,  and  the  remaining  30%  for  validation.  The  available  density  data  were  not  evenly  distributed  in 
time;  therefore,  for  implementing  the  neural  network,  a  linear  interpolation  was  applied  to  make  sure  that 
there  was  a  constant  difference  in  time  between  consecutive  samples  in  the  data.  The  values  of  the  density 
are  in  the  order  of  magnitude  of  10  12  kg/m3  for  day  140  of  2002.  This  results  in  numerical  problems 
during  the  training  of  the  NNs.  To  address  this  issue,  the  natural  logarithm  of  the  density  values  was  used 
for  the  NNs  instead  of  the  density  values  themselves.  Another  advantage  of  using  the  natural  logarithm, 
shown  in  [24],  is  that  it  often  stabilizes  the  variance  of  the  series,  which  allows  for  better  modeling  of  the 
time  series. 

Several  different  periods  of  interest  for  training,  validation,  and  testing  the  NNs  were  identified. 
Stastny  et  al.  [8]  chose  two  representative  days  for  low  and  high  geomagnetic  activities  for  testing  his 
linear  model  and  for  comparing  it  to  JB2006  and  HASDM.  The  first  of  these  days  was  day  141  of  2002; 
during  this  day  there  was  very  low  geomagnetic  activity  (Dst=-16,  ap=10  and  F|07=  190.4).  The  second 
day  used  by  Stastny  et  al.  [8]  was  day  276  of  2001.  During  this  day  there  was  a  moderate  geomagnetic 


storm,  and  as  a  result  there  was  a  higher  geomagnetic  activity  (Dst  =-107,  ap=69  and  F107=191.8)  For 
obtaining  the  linear  model,  Stastny  et  al.  [8]  used  the  data  from  day  140  (Dst  =-12,  ap=10  and 

Fio.7=175.4).  This  same  data  set  was  used  to  train,  validate  and  test  the  NNs.  These  data  sets  included 
77=1080  data  points  for  each  day  with  a  sampling  rate  of  80  sec. 

To  study  the  long-term  performance  of  the  NNs,  it  was  decided  to  test  them  over  one -year  intervals. 
Out  of  the  years  that  CHAMP  was  collecting  data,  years  2003  and  2007  were  of  special  from  the  point  of 
view  of  space  weather  and  therefore  were  selected  for  testing.  During  2003  (Dst  =-22,  ap=21.8  and 
Fio.7=128)  the  geomagnetic  activity  was  the  highest  of  that  solar  cycle  (as  explained  in  [25]).  In  contrast, 
during  2007  (Dst  =-8,  ap=7.5  and  Fi0.7=73)  the  solar  cycle  went  through  a  period  of  very  low  activity 
(solar  minimum,  see  [26])  and  therefore  the  solar  and  geomagnetic  activities  were  very  low.  The  training 
and  validation  sets  consisted  of  data  from  years  previous  to  the  testing  years.  Specifically,  data  from  2002 
was  used  to  train  the  NN  for  2003  prediction,  and  data  from  2006  was  used  to  train  the  NN  for  2007 
prediction.  For  these  long-term  experiments  a  sampling  rate  of  120  sec  was  used.  Figure  3  shows  the  daily 
averaged  values  for  Dst,  ap,  and  Fi0.7  during  2003  and  2007. 
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Figure  3.  Dst,  ap  and  F10.7  indices  averaged  daily  for  years  2003  (left)  and  2007  (right). 

2.2.3  Solar  and  geomagnetic  indices 


By  including  additional  inputs  other  than  the  present  density  values,  the  performance  of  a  NN  as  a 
predictor  improves,  provided  that  the  output  of  the  NN  is  a  function  of  these  inputs.  Because  the  density 
is  driven  by  the  solar  and  geomagnetic  activities,  one  proxy  for  each  of  these  was  selected  as  additional 
inputs.  Out  of  the  many  possible  choices  (Sjo,  M10,  Mg  II)  F10  7  was  assumed  to  be  a  suitable  proxy  for  the 


solar  activity  affecting  the  density.  This  is  a  valid  assumption  since  as  pointed  out  by  [27]  there  has  not 
yet  been  any  proof  of  one  index  been  clearly  superior  to  any  other  for  satellite  operations.  Dst  was 
assumed  to  be  a  valid  proxy  for  geomagnetic  activity.  The  ap  index  was  not  used  since,  as  Figure  3 
shows,  Dst  and  ap  are  closely  related.  Furthermore  as  pointed  out  in  [28],  replacing  ap  for  Dst  reduced 
density  errors  especially  during  geomagnetic  storms  for  Jacchia  70,  NRLMSIS  and  JB  2008. 

The  indices  were  averaged  hourly  and  were  included  in  the  corresponding  training,  validation,  and 
testing  sets.  No  interpolation  was  necessary  for  the  indices  since  their  sampling  rate  (hourly  averaged 
values)  was  much  larger  than  the  sampling  rate  of  the  density  (80  and  120  seconds)  and  there  were  no 
gaps  in  the  data.  For  the  one  orbit  prediction  case  at  a  sampling  rate  of  80  seconds,  68  samples  per 
window  are  used;  for  the  eight  orbits  case  at  a  sampling  rate  of  120  seconds,  360  samples  per  window  are 
used;  and  for  the  32  orbits  case  at  a  sampling  rate  of  120  seconds,  1440  samples  per  window  are  used.  As 
with  the  density  values,  during  operation  the  NNs  only  have  access  to  present  values  of  the  indices.  The 
values  for  the  Dst  and  Fi0  7  indices  used  in  this  work  were  obtained  from  [29]. 


3.  Numerical  Results 


3.1  Maneuver  simulations 

The  three  different  Lyapunov  controllers,  developed  in  the  previous  section  were  tested  and  evaluated 
using  numerical  simulations.  The  guidance  and  control  formulations  have  been  programmed  in 
MATLAB.  These  algorithms  interact  with  STK  via  STK  Connect.  STK’s  High-Precision  Orbit 
Propagator  (HPOP)  was  used  for  modeling  the  mechanics  of  the  maneuver,  including  J2  perturbations, 
solar  pressure  radiation  and  variable  atmospheric  drag  using  the  empirical  NRLMSISE-00  model.  The 
linear  reference  model  and  guidance  trajectories  are  modeled  using  Simulink  in  the  initialization  part  of 
the  simulation  (prior  to  running  STK).  The  simulation  architecture  can  be  seen  in  Figure  4.  For  all  of  the 
simulations,  the  non-adaptive  version  (with  constant  P  and  AJ)  was  compared  with  the  adaptive  and 
optimized  adaptive  Lyapunov  controllers  (with  variable  P  and  AJ).  For  the  re -phase  and  rendezvous 
maneuvers,  the  simulations  were  stopped  when  the  spacecraft  were  within  10  m  of  the  desired  final  state. 
For  the  fly-around  maneuver  the  simulation  was  stopped  after  2.5  orbital  periods  after  the  guidance 
reached  the  final  equilibrium  orbit. 


SIMULINK 

Figure  4.  Simulation  Diagram. 

The  three  controllers  can  be  implemented  in  the  following  configurations: 

1)  The  controller  forces  the  nonlinear  system  to  go  to  a  desired  final  state  (regulation). 

2)  The  controller  is  used  to  force  the  nonlinear  system  to  directly  track  a  generated  guidance 
trajectory  (no  linear  dynamics  is  defined). 

3)  The  controller  forces  the  nonlinear  system  to  track  the  trajectory  of  the  reference  model 
which  is  tracking  the  generated  guidance  trajectory. 


The  initial  orbital  elements  of  the  target  (center  of  the  LVLH  frame)  and  other  parameters  for  the 
numerical  simulations  are  shown  in  Table  1.  The  target  and  chaser  are  assumed  to  be  identical  spacecraft, 
therefore  they  have  the  same  ballistic  coefficient  and  mass,  for  each  surface  configurations.  The  initial 
relative  position  and  velocity  of  the  chaser  in  the  LVLH  frame  are  shown  in  Table  2.  For  all  simulations 
the  initial  Q  matrix  was  the  identity  matrix  times  102. 

Table  1.  Spacecraft  orbital  Parameters  and  geometry. 


Parameter 

Value 

Target’s  inclination  (deg) 

98 

Target’s  semi-major  axis  (km) 

6778 

Target’s  right  ascension  of  the  ascending  node  (deg) 

262 

Target’s  argument  of  perigee  (deg) 

30 

Target’s  true  anomaly  (deg) 

25 

Target’s  eccentricity 

0 

Maximum  Ballistic  Coefficient,  surfaces  deployed  (nr/kg) 

0.625 

Minimum  Ballistic  Coefficient,  surfaces  retracted  (nr/kg) 

0.075 

Mass  (kg) 

10 

Table  2.  Initial  conditions  in  the  LVLH  frame. 


Parameter 

Rendezvous 

Fly-around 

Re-phase 

x  (km) 

-1 

0 

0 

y  (km) 

-2 

-4.25 

-1.9 

X  (km/sec) 

4.8E-07 

0 

0 

y  (km/sec) 

1.70E-04 

0 

0 

These  initial  conditions  where  generated  by  selecting  the  orbits  of  the  target  and  chaser  and  then 
calculating  initial  state  vector  (position  and  velocity  of  the  chaser  relative  to  the  target). 

3.1.1  Re-phase  maneuver 

For  this  simulation  in  STK,  the  initial  difference  in  the  y  was  -1.9km  and  the  desired  final  difference 
was  three  km.  The  three  controllers  (Lyapunov,  adaptive  Lyapunov,  and  optimized  adaptive  Lyapunov 
controllers)  were  used  in  simulations  for  the  re -phase  maneuver.  The  controllers  are  used  to  regulate  the 
error  between  the  simulated  relative  positions  and  velocities  and  the  desired  final  relative  positions  and 
velocities.  For  this  maneuver  the  RLqr  value  used  to  obtain  the  initial  Aj  was  10ls.  The  trajectories  in  the 
LVLH  are  compared  in  Figure  5.  Figure  6  shows  the  control  signals  for  the  three  controllers.  Tracking 
error  plots  are  shown  in  Figure  7. 
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Figure  5.  Re-phase  maneuver  trajectories  in  the  x-y  plane:  (left)  complete  maneuver  and  (right)  final  stages  of  the  maneuver. 
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Figure  6.  Re-phase  maneuver  control  signals:  (top)  non-adaptive  Lyapunov  controller,  (middle)  adaptive  Lyapunov  controller  and 
(bottom)  adaptive  optimized  Lyapunov  controller. 
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Figure  7.  Tracking  error  over  the  entire  re-phase  maneuver:  (top)  x  error  and  (bottom)  y  error. 

3.1.2  Fly-around  maneuver 

Again  the  three  controllers  (Lyapunov,  adaptive  Lyapunov,  and  optimized  adaptive  Lyapunov 
controllers)  have  been  used  in  simulations  for  the  fly-around  maneuver.  In  this  case  the  controllers  force 
the  nonlinear  dynamics  to  follow  the  desired  fly-around  guidance  and  not  just  to  converge  to  a  final  state 
(second  configuration).  Moreover,  the  simulations  are  stopped  2.5  orbital  periods  after  the  guidance 
reaches  the  final  equilibrium  orbit.  For  this  maneuver  the  RLqr  value  was  1 .6*  1 018.  The  guidance  selected 
for  the  fly-around  maneuver  is  based  on  Clohessy -Wiltshire  equations  and  was  obtained  from  [30].  The 
trajectories  in  the  LVLH  are  compared  in  Figure  8,  while  Figure  9  shows  the  control  signals  for  the  three 
controllers.  Tracking  error  plots  are  shown  Figure  10. 
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Figure  8.  Fly-around  trajectories  in  the  x-y  plane:  (left)  complete  maneuver  and  (right)  final  stages  of  the  maneuver. 
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Figure  9.  Fly-around  maneuver  control  signals:  (top)  non-adaptive  Lyapunov  controller,  (middle)  adaptive  Lyapunov  controller 
and  (bottom)  adaptive  optimized  Lyapunov  controller. 
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Figure  10.  Tracking  error  over  the  entire  fly-around  maneuver:  (top)  x  error  and  (bottom)  y  error. 

3.1.3  Rendezvous  maneuver 

Of  the  three  maneuvers  simulated,  the  rendezvous  maneuver  is  the  most  challenging  using  DD  since  it 
requires  reducing  an  initial  difference  in  the  x  direction  as  well  as  in  the  y  direction.  Since  the  system  is 
underactuated,  this  requires  introducing  large  errors  in  the  y  direction.  For  this  maneuver  the  controllers 
have  been  implemented  using  the  reference  model  tracking  (third  configuration). 

A  guidance  trajectory  based  on  Clohessy-Wiltshire  equations  has  been  selected  for  the  rendezvous 
maneuvers.  This  guidance  trajectory  was  made  following  the  method  described  in  [3]  which  uses  a 
constant  value  for  the  density.  In  this  simulation  the  controller  forces  the  nonlinear  system  to  track  the 
trajectory  of  the  reference  model  which  is  tracking  the  analytically  generated  guidance  trajectory  (third 
configuration).  For  this  maneuver  the  RLqr  value  was  1.6*1017.  The  trajectories  in  the  LVLH  are 
compared  in  Figures  11,  12  and  13.  Figure  14  shows  the  control  signals  for  the  three  controllers.  Tracking 
error  plots  are  shown  Figure  15. 
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11.  Rendezvous  trajectories  for  the  non-adaptive  controller  in  the  x-y  plane:  (left)  complete  maneuver  and  (right)  final 
of  the  maneuver. 
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Figure  12.  Rendezvous  trajectories  for  the  adaptive  controller  in  the  x-y  plane:  (left)  complete  maneuver  and  (right)  final  stages 
of  the  maneuver. 
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Figure  13.  Rendezvous  trajectories  for  the  optimized  adaptive  controller  in  the  x-y  plane:  (left)  complete  maneuver  and  (right) 
final  stages  of  the  maneuver. 
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Figure  14.  Rendezvous  maneuver  control  signals:  (top)  non-adaptive  Lyapunov  controller,  (middle)  adaptive  Lyapunov  controller 
and  (bottom)  adaptive  optimized  Lyapunov  controller. 
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Figure  15.  Tracking  error  over  the  entire  rendezvous  maneuver:  (top)  x  error  and  (bottom)  y  error. 

3.1.4  Controller  performance  assessment 

The  metrics  used  to  evaluate  the  performance  of  the  controllers  are  the  number  of  switches  in  the 
control  (control  effort),  the  duration  of  the  maneuver,  the  means  for  the  critical  and  actual  value  of  the  DD 
acceleration,  and  the  difference  between  these  two  values  (control  margin).  The  performance  metrics  for 
the  simulations  are  shown  in  Table  3. 


Table  3.  Performance  metrics  for  all  maneuvers 


Metric 

Re-phase 

(Regulation) 

Fly-around 

(Tracking 

Trajectory) 

Rendezvous  (Tracking 
Dynamics) 

Control 

changes 

124 

48 

79 

Time  (hr) 

30.7 

13.2 

43.6 

Non- 

Adaptive 

Drag  Mean 
Critical  value 
(m/s2) 

-2.93E-07 

-5.91E-06 

-1.43E-04 

Mean  Actual 
Drag(m/s2) 

9.61E-06 

3.38E-05 

3.42E-05 

Margin(m/s2) 

9.90E-06 

3.97E-05 

1.77E-04 

Control 

changes 

107 

45 

37 

Time  (hr) 

27.4 

13.2 

35.4 

Adaptive 

Drag  Mean 
Critical  value 
(m/s2) 

-6.54E-07 

-6.36E-06 

-1.81E-04 

Mean  Actual 
Drag(m/s2) 

9.52E-06 

3.38E-05 

3.49E-05 

Margin(m/s2) 

1.02E-05 

4.01E-05 

2.16E-04 

Control 

changes 

73 

46 

35 

Time  (hr) 

19.8 

13.2 

35.1 

Optimized 

Adaptive 

Drag  Mean 
Critical  value 
(m/s2) 

-2.05E-06 

-6.12E-06 

-1.82E-04 

Mean  Actual 
Drag(m/s2) 

9.44E-06 

3.38E-05 

3.46E-05 

Margin(m/s2) 

1.15E-05 

3.99E-05 

2.17E-04 

As  shown  in  Figures  5,  7,  11,  12,  13  and  15  for  the  re -phase  and  rendezvous  maneuvers,  all  three 
controllers  forced  the  chaser  and  target  spacecraft  to  be  within  10m  of  their  desired  final  state.  As  seen  in 
see  Figure  8  for  the  fly-around  maneuver  the  two  adaptive  controllers  forced  the  system  to  reach  an  orbit 
close  to  the  desired  equilibrium  final  orbit  (within  ~50m);  however,  the  non-adaptive  reached  an  orbit 
with  less  precision  than  the  other  controllers  (within  ~300m).  Furthermore,  as  can  be  observed  in  the 
tracking  error  plots  (Figures  7,  10  and  15)  both  adaptive  controllers  produce  smaller  oscillations  on  the 
error,  especially  toward  the  end  on  the  maneuvers,  this  produces  smoother  trajectories;  even  though,  the 
control  force  is  varying  (since  the  density  is  varying). 

The  use  of  the  adaptive  controller  reduced  the  control  effort  required  to  perform  all  three  maneuvers 
with  improvements  relative  to  the  non-adaptive  of  13.7%,  6.3%,  and  53.2%  for  the  re -phase,  fly-around 
and  rendezvous  maneuvers  respectively.  Similarly,  the  use  of  the  optimized  adaptive  controller  reduced 


the  control  effort  by  41.1%,  4.2%,  and  55.7%  in  comparison  with  the  non-adaptive  for  the  re -phase,  fly- 
around  and  rendezvous  maneuvers  respectively. 

Likewise,  the  use  of  adaptive  controller  also  reduced  the  duration  of  the  maneuver  with  improvements 
of  10.9%,  and  18.8%  for  the  re -phase  and  rendezvous  maneuvers  respectively.  Similarly,  the  use  of  the 
optimized  adaptive  controller  reduced  the  duration  of  the  maneuver  with  improvements  of  35.7%,  and 
19.6%  over  the  non-adaptive  for  the  re -phase  and  rendezvous  maneuvers  respectively.  There  were  no 
improvements  in  duration  for  the  fly-around,  since  the  simulations  were  not  stopped  when  within  10m  of 
the  desired  final  position,  but  after  2.5  orbital  periods  after  reaching  the  desired  equilibrium  relative  orbit. 

Furthermore,  the  adaptive  controller  also  was  able  to  increase  the  control  margin  by  2.9%,  1%,  and 
18.1%  over  the  non-adaptive  for  the  re -phase,  fly-around,  and  rendezvous  maneuvers  respectively. 
Similarly,  the  optimized  adaptive  controller  increased  the  control  margin  over  the  non-adaptive  by  13.9%, 
0.5%,  and  18.4%  over  the  non-adaptive  for  the  re -phase,  fly-around,  and  rendezvous  maneuvers 
respectively. 

Overall  both  adaptive  and  optimized  adaptive  controllers  gave  better  results  than  the  non-adaptive 
controller.  The  optimized  adaptive  gave  better  results  than  the  adaptive  for  the  re -phase  (Regulation)  and 
rendezvous  (Tracking  Dynamics)  maneuvers;  however,  for  the  case  of  the  fly-around  (Tracking 
Trajectory),  the  performance  metrics  of  the  optimized  adaptive  were  slightly  worse  than  the  adaptive. 
Nonetheless,  Figure  10  shows  that  during  the  maneuver  the  errors  were  smaller  for  the  optimized 
adaptive,  so  perhaps  this  came  at  a  cost  of  higher  slightly  higher  actuation  (one  more  switch  in  the  drag 
surfaces  configuration). 

3.2  Density  Predictors 

The  training,  validation,  and  testing  of  the  NNs  was  done  in  MATLAB  using  the  Neural  Network 
Toolbox.  As  a  benchmark  for  all  the  tests,  a  model  using  the  persistence  method  was  used.  The 
persistence  method  is  a  very  simple  technique  for  forecasting  in  which  the  prediction  is  equal  to  the  input. 
In  other  words,  a  predictor  forecasting  the  density  using  the  persistence  method  will  predict  the  density  in 
the  future  to  be  the  same  as  it  is  in  the  present,  so  if  the  prediction  window  is  one  orbit,  then  each 
predicted  orbit  is  equal  to  the  previous  measured  orbit. 

To  assess  the  performance  of  the  different  models,  different  metrics  were  used:  the  MSE  (shown  in 
Equation  (15)),  the  mean  of  the  ratio  between  the  target  and  the  outputs,  its  standard  deviation  and  the 
Pearson  correlation  coefficient  of  the  targets  to  the  model  outputs  ( Rp ). 


3.2.1  Different  number  of  neurons,  delays  and  sampling  rates 

Given  the  selected  structure  of  the  NN  predictors,  the  appropriate  number  of  neurons,  delays  in  the 
hidden  layer,  and  the  data  sampling  rate  for  the  localized  density  forecasting  problem  were  found.  This 
was  accomplished  empirically  by  trying  different  combinations.  All  the  tests  performed  for  this  purpose 
were  run  on  days  141  of  2002  and  276  of  2001,  with  the  training  and  validation  sets  being  day  140  of 
2002.  As  mentioned  before,  these  days  cover  high  and  low  geomagnetic  activity  and  were  used  also  by 
Stastny  et  al.  [8]  to  test  his  linear  model.  To  find  the  appropriate  number  of  neurons  and  delays  in  the 
hidden  layer  and  the  sampling  rate  the  NN  shown  in  Figure  16  was  used. 
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Figure  16.  Neural  network  diagram  where  N  and  D  are  the  numbers  of  neurons  and  delays  in  the  hidden  layer. 

To  find  the  appropriate  number  of  neurons  in  the  hidden  layer  several  tests  were  performed  in  the 
sampling  rate  and  the  number  of  delays  were  fixed  and  the  number  of  neurons  was  varied.  Afterwards,  the 
same  procedure  was  followed  for  the  number  of  delays  and  the  sampling  rate.  The  results  for  these  tests 
are  presented  in  [15].  The  best  results  were  obtained  with  one  neuron  in  the  hidden  layer  and  enough 
delays  to  store  from  lA  to  one  prediction  window,  while  there  was  no  significant  difference  in  the 
performance  when  the  sampling  rate  was  below  180  seconds. 

3.2.2  Predicting  one  orbit  into  the  future  on  days  241  of  2002  and  276  of  2001 

Once  the  appropriate  structure  of  the  NN  was  found  several  different  NNs  were  tested  again  on  days 
141  of  2002  and  276  of  2001.  This  was  done  to  evaluate  the  improvements  in  performance  by  increasing 
the  size  of  the  training  and  validation  sets  from  one  day  to  a  year  and  also  by  using  the  solar  and 
geomagnetic  indices  (Dst  and  F10.7)  as  additional  inputs.  The  one  year  training  and  validation  data  set  used 
for  testing  the  networks  on  day  141  of  2002  contained  the  data  from  the  365  preceding  days  (day  140  of 
2001  to  day  14  of  2002).  The  one  year  training  and  validation  data  set  used  for  testing  the  networks  on 
day  276  of  2001  contained  the  data  from  year  2002  (day  1  of  2002  to  day  365  of  2002),  since  the 
CHAMP  data  did  not  went  back  a  year  before  day  276  of  2001 .  Even  though  this  means  that  the  NN  used 
was  trained  on  data  corresponding  to  the  future  of  day  276  of  2001,  the  training  data  and  validation  data 
set  is  still  different  to  the  testing  set  which  makes  the  test  valid  (of  course  for  practical  implementation  of 
the  NNs  the  training  and  validation  set  would  always  be  past  and  therefore  available  values).  A  sampling 


rate  of  80  sec  was  used  since  it  is  the  same  used  by  Stastny  et  al.  [8].  The  results  of  the  tests  for  days  141 
of  2002  and  276  of  2001  are  summarized  in  Table  4.  As  a  benchmark  the  performance  of  a  persistence 
model  are  also  included.  The  persistence  method  is  a  very  simple  technique  for  forecasting  in  which  the 
prediction  is  equal  to  the  input. 

Table  4.  Results  for  predicting  one  orbit  into  the  future. 


Testing 
Data  Set 

Model 

Configuration 

MSE 

Mean 

target/output 

Stdev 

target/output 

NN,  preceding 

365  days  for 
Training 

0.0108 

0.9843 

0.9998 

0.0039 

Day  141 
2002 
CHAMP 

NN,  preceding 

365  days  for 
Training,  Dst  and 
F|0.7 

0.0108 

0.9842 

0.9998 

0.0039 

NN,  day  140  of 
data  for  Training 

0.0156 

0.9774 

1.0008 

0.0046 

Persistence  Model 

0.0234 

0.9685 

1.0003 

0.0058 

Linear  model* 

N/A 

N/A 

1.0058 

0.0822 

HASDM* 

N/A 

N/A 

0.8662 

0.1204 

JB2006* 

N/A 

N/A 

0.8564 

0.095 

NN,  year  2002  of 
data  for  Training 

0.0229 

0.9086 

0.9999 

0.0058 

NN.  Year  2002  of 
data  for  Training, 
Dst  and  Fi0  7 

0.0225 

0.9099 

0.9999 

0.0058 

Day  276 
2001 
CHAMP 

NN,  day  140  Year 
2002  of  data  for 
Training 

0.0229 

0.9106 

0.9995 

0.0058 

Persistence  Model 

0.0328 

0.8718 

1.0001 

0.007 

Linear  model* 

N/A 

N/A 

1.0094 

0.0822 

HASDM* 

N/A 

N/A 

0.8415 

0.1344 

JB2006* 

N/A 

N/A 

0.6471 

0.1355 

*  Obtained  from  [8] 

For  the  2001  scenario,  training  data  from  2002  is  used,  as  Stastny  et  al.  did  in  their  work.  Training 
with  future  values  and  "predicting"  past  values  is  valid  from  the  point  of  view  of  neural  network,  since  the 
training/validation  and  testing  data  sets  are  still  different.  The  results  in  Table  4  indicate  that  the  global 
models  (HASDM  and  JB  2006  results  obtained  from  Stastny  et  al.  [8])  have  large  biases  in  their  results 
for  the  test  days.  This  causes  their  performance  to  be  worse  than  the  performance  of  all  the  other  models 
including  the  persistence  model.  The  NN  predictors  give  significantly  better  results  than  the  linear  model 
from  Stastny  et  al.  [8]),  the  global  models,  and  the  persistence  model.  For  day  141  of  2002,  by  increasing 


the  size  of  the  training  and  validation  sets,  the  performance  of  the  NNs  increases;  however,  for  day  276  of 
2001  there  is  not  a  significant  improvement  by  increasing  the  size  of  the  training  and  validation  sets  nor 
by  including  the  solar  and  geomagnetic  indices.  The  addition  of  the  indices  does  not  benefit  the  NNs 
because  the  number  of  delays  (17  which  corresponds  to  lA  of  the  prediction  window)  cannot  capture  more 
than  one  value  in  time  of  the  indices  since,  the  indices  are  averaged  hourly.  This  might  be  solved  by 
increasing  the  number  of  delays.  An  alternative  solution  would  be  retaining  the  same  number  of  delays, 
but  space  them  non-uniformly  in  time. 

For  day  141  of  2002,  using  data  from  one  year  to  train  and  validate  the  NN  provided  the  best  results. 
The  actual  output  of  this  NN  and  the  targets  are  shown  in  Figure  17  along  with  the  prediction  error.  For 
day  276  of  2001,  the  NN  that  uses  the  additional  inputs  (Dst  and  F10.7)  and  that  was  trained  and  validated 
using  the  data  from  one  year  yielded  the  best  results.  The  actual  output  of  this  neural  network,  the  targets, 
and  the  prediction  error  are  shown  in  Figure  18. 
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Figure  17.  Neural  network  response  for  best  case  with  a  prediction  window  of  one  orbit  over  day  141  of  2002. 
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Figure  18.  Neural  network  response  for  best  case  with  a  prediction  window  of  one  orbit  over  day  276  of  2001. 


3.2.3  Predicting  eight  and  32  orbits  into  the  future 

For  most  applications  of  the  NN  density  predictors,  longer  prediction  windows  are  desired.  For  this 
reason  additional  NN  predictors  were  trained,  validated,  and  tested  for  predicting  eight  and  32  orbits  into 
the  future  (approximately  half  a  day  and  two  days).  For  these  results,  the  NNs  were  tested  on  years  2003 
and  2007,  in  order  to  evaluate  their  performance  over  much  wider  data  sets  including  periods  of  low  and 
high  solar  and  geomagnetic  activities.  As  before,  the  use  of  additional  inputs  (Dst  and  Fi0  7)  was  studied 
along  with  the  use  of  different  numbers  of  delays.  Since  having  different  sampling  rates,  as  long  as  they 
are  below  180  sec,  does  not  affect  significantly  the  NN  performance,  a  sampling  rate  of  120  sec  was  used 
in  order  to  reduce  time  for  training  and  validation  the  NNs.  The  results  are  summarized  in  Tables  5  and  6 
for  the  prediction  windows  of  eight  and  32  orbits  respectively. 


Table  5.  Results  for  predicting  eight  orbits  into  the  future. 


Testing 
Data  Set 

Model 

Configuration 

MSE 

Mean 

target/output 

Stdev 

target/output 

ANN,  90  delays  (2 
orbits),  1  year  2002 
of  data  for  Training 

0.0433 

0.8971 

1.0007 

0.0078 

CHAMP 

2003 

ANN,  90  delays  (2 
orbits),  1  year  2002 
of  data  and  Dst  and 
F10.7  for  Training 

0.0429 

0.8976 

1 

0.0078 

ANN,  360  delays  (8 
orbits),  1  year  2002 
of  data  and  Dst  and 
F10.7  for  Training 

0.0401 

0.9044 

0.9999 

0.0075 

Persistence  Model 

0.2614 

0.4037 

1.0002 

0.0192 

ANN,  90  delays  (2 
orbits),  1  year  2006 
of  data  for  Training 

0.0417 

0.9093 

1.0002 

0.0075 

CHAMP 

2007 

ANN,  90  delays  (2 
orbits),  1  year  2006 
of  data  and  Dst  and 
Fjo  7  for  Training 

0.0407 

0.9114 

1 

0.0074 

ANN,  360  delays  (8 
orbits),  1  year  2006 
of  data  and  Dst  and 
F10.7  for  Training 

0.0403 

0.9122 

1 

0.0074 

Persistence  Model 

0.1902 

0.6031 

1.0001 

0.016 

The  best  case  included  in  Table  5  for  both  years  2003  and  2007  were  those  obtained  with  the  NN  that 
included  the  additional  inputs  and  that  had  360  delays  (one  prediction  window).  Figures  19  and  20  show 
the  MSE  over  the  entire  years  2003  and  2007  for  the  best  cases  along  with  the  Dst  and  F10.7  averaged 
daily. 


MSE 


MSE  over  one  year 


Figure  19.  MSE  for  best  case  with  a  prediction  window  of  eight  orbits  and  normalized  indices  over  year  2003. 
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Figure  20.  MSE  for  best  case  with  a  prediction  window  of  eight  orbits  and  normalized  indices  over  year  2007. 
Table  6.  Results  for  predicting  32  orbits  into  the  future. 


Testing 
Data  Set 

Model 

Configuration 

MSE 

Mean 

target/output 

Stdev 

target/output 

ANN,  360  delays  (8 
orbits),  1  year  2002  of 
data  for  Training 

0.0917 

0.7702 

1.0013 

0.0113 

CHAMP 

2003 

ANN,  360  delays  (8 
orbits),  1  year  2002  of 
data  and  Dst  and  F10,7 
for  Training 

0.0895 

0.774 

1 

0.0112 

Persistence  Model 

0.1813 

0.5874 

1.0001 

0.016 

ANN,  360  delays  (8 
orbits),  1  year  2006  of 
data  for  Training 

0.1564 

0.6058 

1.0009 

0.0145 

CHAMP 

2007 

ANN,  360  delays  (8 
orbits),  1  year  2006  of 
data  and  Dst  and  FI0.7 
for  Training 

0.1515 

0.6215 

1.0002 

0.0143 

Persistence  Model 

0.603 

-0.2582 

1.0005 

0.0286 

As  indicated  in  Table  6  for  the  32  orbit  predictions  the  NNs  yielded  a  better  performance  over  the  high 
activity  period  (2003)  than  the  low  activity  one  (2007).  This  was  not  observed  in  any  of  the  other  tests 
performed  where  the  performance  for  both  periods  was  almost  the  same  or  better  for  periods  of  low 
activity  (see  Tables  4  5).  This  indicates  that  for  longer  prediction  windows  over  periods  of  low  activity 
other  unknown  factors  affect  the  density  behavior,  that  are  not  well  represented  by  the  data  used  by  the 
NNs  (current  value  of  the  density,  Dst  and  Fi0.7  indices).  This  is  here  considered  a  topic  for  further 
investigation  and  beyond  the  scope  of  this  work.  Further  investigation  may  lead  to  the  discovery  of 
unknown  effects  during  periods  of  low  activity. 

The  best  case  shown  in  Table  6  for  both  years  2003  and  2007  were  those  obtained  with  the  NN  that 
used  the  Dst  and  Fi0  7  indices.  The  MSE  over  the  entire  years  2003  and  2007  for  the  best  cases  included  in 
Table  6  (32  orbits  prediction),  along  with  the  Dst  and  F10.7  averaged  daily  are  shown  in  Figures  21  and  22 
for  years  2003  and  2007  respectively. 


MSE  over  one  year 


Figure  21.  MSE  for  best  case  with  a  prediction  window  of  32  orbits  and  normalized  indices  over  year  2003. 
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One  very  interesting  feature  can  be  observed  in  Figures  19  and  21.  The  peaks  in  the  MSE  correspond 
to  peaks  in  the  Dst  index.  This  indicates  that  the  NN  predictors  will  have  larger  errors  during  geomagnetic 
storms.  This  is  further  confirmed  by  the  results  from  the  test  performed  for  predicting  one  orbit  into  the 
future,  in  which  the  performance  for  the  predictions  on  day  141  of  2002  were  always  better  than  the 
performance  on  day  276  of  2001,  during  which  there  was  a  geomagnetic  storm. 

From  a  computational  point  of  view  the  training  and  validation  of  the  NN  predictors  can  be  costly 
especially  when  dealing  with  large  prediction  windows  (32  orbits  into  the  future).  For  this  reason,  for 
future  implementation  of  this  work  the  training  and  validation  processes  are  not  recommended  to  be  done 
onboard.  Rather,  it  is  proposed  that  density  values  obtained  onboard  via  accelerometers  be  sent 
periodically  to  the  mission  team  on  the  ground  for  training  and  validation.  Once  the  weights  and  biases  of 
the  trained  NN  are  obtained,  they  can  be  uplinked  to  the  onboard  computers  and  then  the  predictor  can  be 
used  for  onboard  orbit  propagation.  The  NN  could  then  be  re -trained  on  the  ground  as  necessary. 


4.  Conclusions 


This  work  developed  methods  for  controlling  the  relative  motion  of  a  target  and  chaser  spacecraft  in 
the  orbital  plane  using  a  differential  in  their  drag  accelerations,  to  enable  propellant-less  coplanar  relative 
maneuvers  in  LEO.  By  varying  the  DD  between  the  chaser,  and  target  spacecraft  their  relative  motion  in 
the  orbital  plane  is  controlled.  These  variations  are  induced  by  changing  the  crosswind  area  of  the 
spacecraft.  The  nature  of  this  variation  is  assumed  to  be  of  bang-off -bang  nature  with  only  three  possible 
values:  maximum  differential  acceleration,  minimum  differential  acceleration,  and  zero  differential 
acceleration. 

Lyapunov  principles  are  used  to  develop  a  control  law  for  the  activation  of  the  drag  surfaces,  thus 
allowing  for  any  co-planar  propellant -less  spacecraft  relative  maneuvering,  physically  realizable  with  the 
small  relative  accelerations  created  by  the  difference  in  the  drag  forces  acting  on  the  spacecraft.  The 
Lyapunov  controller  can  be  implemented  by  forcing  the  nonlinear  dynamics  to  track  a  trajectory,  the 
dynamics  of  a  reference  model,  or  to  regulate  to  a  desired  final  state.  This  allows  for  the  implementation 
of  the  controller  in  maneuvers  in  which  a  specific  path  is  desired,  consequently,  opening  the  possibilities 
for  any  other  maneuver,  provided  that  they  are  confined  to  the  orbital  plane  and  that  they  are  realizable 
using  the  small  acceleration  created  by  the  drag.  The  analytical  nature  of  the  methodology  holds  the 
promise  for  future  onboard  implementation  on  real  spacecraft. 

An  analytical  expression  the  magnitude  of  the  DD  that  ensures  stability  (critical  value)  was  developed. 
Analytical  expressions  for  the  partial  derivatives  of  the  DD  critical  value  ensuring  Lyapunov  stability  in 
terms  of  matrices  Q  and  Aa  (chosen  by  the  control  designer,  i.e.,  independent  variables),  were  also 
developed.  Two  adaptations  to  the  Lyapunov  control  law  were  developed.  These  adaptations  use  the 
partial  derivative  of  the  critical  value  in  terms  of  matrix  Ad  to  find  the  element  of  this  matrix  to  which  the 
critical  value  is  the  most  sensitive.  In  the  first  adaptation,  this  element  was  changed  by  switching  its  value 
between  fixed  values  to  adapt  the  control  law.  In  the  second  adaptation  this  element  was  changed  to  be 
the  value  that  minimized  (within  a  fixed  range)  the  analytical  expression  of  the  critical  value.  The 
adaptations  result  in  a  matrix  4 A  that  varies  in  time,  and  hence,  so  does  the  P  matrix  in  the  quadratic 
Lyapunov  function.  Thus  a  new  term,  containing  the  time  derivate  of  the  P  matrix,  appears  in  the 
Lyapunov  function.  An  analytical  expression  for  the  time  derivative  of  the  P  matrix  (assuming  a 
continuous  adaptation)  was  obtained  in  order  to  evaluate  the  impact  of  this  new  term  in  the  stability  of  the 
system. 

STK’s  HPOP  was  used  to  simulate  the  nonlinear  dynamics  of  relative  motion  for  spacecraft.  In 
simulations  the  results  for  the  re -phase,  fly-around,  and  the  rendezvous  maneuvers  indicate  that  the 
implementation  of  both  the  adaptive  and  optimized  adaptive  Lyapunov  controllers  allows  for  smoother 
maneuvers  with  less  duration,  less  actuation,  and  greater  control  margin  for  the  three  different  controller 


configurations  studied  with  respect  to  the  non-adaptive  controller.  Furthermore,  the  optimized  adaptive 
offers  slightly  better  performance  than  the  adaptive  in  terms  of  maneuver  duration  and  control  effort.  The 
accuracy  achieved  with  the  proposed  controllers  is  unprecedented  for  DD  based  relative  maneuvering.  For 
the  rendezvous  and  re -phase  maneuvers  all  controllers  had  the  ability  to  take  the  spacecraft  to  within  10 
meters  of  their  desired  final  states,  using  only  the  differential  in  drag.  Furthermore,  for  the  fly  around 
maneuver  both  adaptive  controllers  reached  relative  orbits  very  close  to  the  desired  final  equilibrium 
relative  orbit. 

The  NN  predictors  provided  significantly  better  results  than  a  linear  model,  and  the  global  models 
HASDM  and  JB2006  for  predicting  the  value  of  the  density  one  orbit  into  the  future  for  periods  of  high 
and  low  geomagnetic  activity.  The  NN  predictors  were  also  tested  for  predicting  eight  and  32  orbits  into 
the  future  (about  half  a  day  and  two  days  at  CffAMP’s  orbit).  For  these  tests  the  performance  of  the  NN 
predictors  was  evaluated  over  the  years  2003  and  2007,  which  cover  the  periods  with  the  high  and  low 
solar  and  geomagnetic  activities.  The  performance  of  the  NN  predictors  decreases  as  the  prediction 
window  increases,  but  even  for  the  32  orbit  case,  the  results  were  satisfactory. 

The  NN  predictors  can  also  use  the  current  value  of  the  Dst  (geomagnetic  activity)  and  F10.7  (solar 
activity)  indices  averaged  hourly  as  additional  inputs  which  resulted  in  an  improvement  of  the 
performance  of  the  NNs,  provided  that  enough  delays  were  included  in  the  hidden  layer  to  store  some  of 
the  behavior  in  time  of  the  indices.  However,  the  number  of  delays  cannot  be  increased  beyond  those 
required  to  store  one  prediction  window  or  the  NNs  will  suffer  from  overfitting  in  terms  of  the  density 
values. 

The  controllers  and  the  NN  predictors  are  computationally  simple  and  can  be  implemented  onboard 
spacecraft.  This  would  allow  for  accurate  relative  motion  control  using  DD  and  precise  onboard  orbit 
propagation  that  can  be  used  for  navigation. 
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